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Abstract 

Gap-junctional coupling is an important way of communication between neurons and other excitable 
cells. Strong electrical coupling synchronizes activity across cell ensembles. Surprisingly, in the pres- 
ence of noise synchronous oscillations generated by an electrically coupled network may differ qualita- 
tively from the oscillations produced by uncoupled individual cells forming the network. A prominent 
example of such behavior is the synchronized bursting in islets of Langerhans formed by pancreatic 
/3— cells, which in isolation are known to exhibit irregular spiking [65 , 64 1. At the heart of this intriguing 
phenomenon lies denoising, a remarkable ability of electrical coupling to diminish the effects of noise 
acting on individual cells. 

In this paper, building on an earlier analysis of denoising in networks of integrate-and-fire neurons 
1501 and our recent study of spontaneous activity in a closely related model of the Locus Coeruleus 
network (551, we derive quantitative estimates characterizing denoising in electrically coupled networks 
of conductance-based models of square wave bursting cells. Our analysis reveals the interplay of the 
intrinsic properties of the individual cells and network topology and their respective contributions to 
this important effect. In particular, we show that networks on graphs with large algebraic connectivity 
(201 or small total effective resistance |3 1 are better equipped for implementing denoising. As a by- 
product of the analysis of denoising, we analytically estimate the rate with which trajectories converge 
to the synchronization subspace and the stability of the latter to random perturbations. These estimates 
reveal the role of the network topology in synchronization. The analysis is complemented by numerical 
simulations of electrically coupled conductance-based networks. Taken together, these results explain 
the mechanisms underlying synchronization and denoising in an important class of biological models. 



1 Introduction 

Cells in the nervous system are organized in complex interconnected networks, which feature a rich variety 
of electrical activity. There is abundant experimental evidence linking spatio-temporal features of the firing 
patterns generated by neuronal networks to various physiological and cognitive processes. Therefore, eluci- 
dating dynamical principles underlying pattern-formation in neuronal networks is an important problem of 
mathematical physiology. 
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Square wave bursting, one of the most common firing patterns, is characterized by alternating periods of 
fast oscillations of the membrane potential and quiescence |[58l[2l|. Typically, cells that generate bursting, 
under different conditions exhibit other firing patterns such as periodic or aperiodic spiking or reside in 
the excitable regime ||5l [531 EH [201 • Transitions between different dynamical regimes in excitable cells 
often signal important physiological or cognitive events such as changes in the rate of hormone secretion or 
neurotransmitter release as in the cases of pancreatic /3— cells |65 1 and midbrain dopamine neurons fT7\ : or 
changes in respiratory rhythm or attentional state as in the cases of neurons in the Pre-Botzinger complex H 
or Locus Coeruleus | 72 1 respectively. Not surprisingly, mathematical models have been used extensively to 
explain the origins of different firing patterns. For single cell models, mechanisms underlying various modes 
of electrical activity have been thoroughly studied |[4D[2ni5Tl|58l[59l[69l. Some of the techniques developed 
for single cell models extend to cover small networks ||46l[60l. However, mathematical analysis of large 
conductance-based networks without special assumptions on network topology is an outstanding problem. 
Furthermore, there is a growing body of experimental and theoretical studies indicating the importance of 
noise in shaping neuronal dynamics |[IOl[l5l[23[3lllM[32|42|4ll|^ presence of noise 

in the network dynamics, its analysis becomes even a more challenging problem. The goal of this paper is 
to elucidate principal factors shaping synchronous activity in large electrically coupled networks of bursting 
capable cells forced by small white noise. 

Many ingredients contribute to the output of neuronal networks. Among them, intrinsic properties of the 
individual cells (local dynamics) and the type and the structure of connections between cells (network topol- 
ogy) are probably the most important ones. Direct electrical coupling through gap-junctions is a common 
way of communication between neurons as well as between cells of the heart, pancreas, and other physiolog- 
ical systems [11 J. The role of electrical coupling in shaping firing patterns generated by neuronal networks 
has been studied using many different techniques: the theory for weakly connected networks Il32ll35l[62l , 
Poincare maps ||29l[6l|38lEl, and Lyapunov functions ll54l , to name a few. In the present study, we consider 
a relatively less studied case of strong electrical coupling |[l2l|54 |, for which we develop two complementary 
approaches based on center-manifold reduction ll37ll and fast-slow decomposition ||Il[56l. Importantly, our 
method covers networks with arbitrary topology, which allows us to study a large class of models and to 
reveal the role of the network topology in shaping network dynamics. 

Under fairly general conditions, strong electrical coupling synchronizes activity across the network ||48l 
l49l . Therefore, one might expect that dynamics of electrically coupled networks of bursting cells will closely 
resemble that of a single cell provided the coupling is strong enough. This is true in general for deterministic 
models. However, in the presence of noise network dynamics, while still synchronous, can be qualitatively 
different from that of a single cell uncoupled from its neighbors. For instance, single cell models, which 
exhibit irregular spiking in isolation (Fig. [T^) can generate very regular synchronous bursting when they are 
coupled electrically (Fig. [TJ)). Likewise, a coupled network exhibiting synchronous spiking for extremely 
long period of time (Fig. [T]:) may be formed from bursting cells (Fig. [T]l). The first scenario illustrated 
by Fig. [T^,b was proposed in (651 [64! to explain why pancreatic /3— cells burst within electrically coupled 
islets of Langerhans, but in isolation exhibit irregular spiking. Numerical experiments and formal analysis 
in [65l|64l show that noise shaping the dynamics of individual /3— cells becomes less effective when these 
cells are coupled electrically. It is certainly intuitive (albeit not obvious mathematically) that in a coupled 
network the effects of uncorrelated stochastic processes acting on individual cells may become weaker due to 
averaging. However, a remarkable property of electrical coupling, which as the analysis below shows should 
not be taken for granted, is that the variability of the coupled system can be fully controlled by varying two 
parameters: the coupling strength and the size of the network. We identify analytical conditions, which 
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Figure 1: Scenario A: Irregular spiking generated by the single models ( |2.1| )-( |2.3| ) shown in (a) is trans- 
formed into bursting as shown in (b) when the cells are coupled electrically. A complementary Scenario B is 
illustrated in (b) and (c). The single cell model generates irregular bursting (c). After the coupling is turned 
on the pattern of firing is effectively switched to spiking (d). 



guarantee that the variability of the coupled systems (the meaning of variability will be explained below) 
goes to zero as the network size and the coupling strength tend to infinity. 

Besides bursting in electrically coupled islets of Langerhans, other phenomena where denoising plays a 
central role, include episodes of phasic firing in the Locus Coeruleus network 11721 and enhanced reliability 
of neural responses in gap-junctionally coupled networks ||50l[68l- In this paper, building on an earlier anal- 
ysis of denoising in networks of integrate-and-fire neurons lISOll and our recent study of spontaneous activity 
in a closely related model of the Locus Coeruleus network |[55]| , we derive quantitative estimates charac- 
terizing denoising in electrically coupled networks of conductance-based models. We find that the results 
obtained for integrate-and-fire models for individual cells do not extend automatically to conductance-based 
models with higher-dimensional state phase. We identify additional features of the local dynamics and cou- 
pling architecture that are needed to guarantee denoising. In particular, our analysis highlights the role of 
the bifurcation structure of the bursting cell model for denoising. It also elucidates the contribution of the 
network topology to this important effect. We show that networks on the graphs with large algebraic con- 
nectivity EOl or small total effective resistance (3] are better equipped for implementing denoising. As a 
by-product of the analysis of denoising, we analytically estimate the rate with which trajectories converge 
to the synchronization subspace and the stability of the latter to random perturbations. Taken together, these 
results explain the mechanisms underlying synchronization and denoising in an important class of biological 
models. 

The organization of the paper is as follows. In the next section, we formulate our assumptions on the 
single cell model ( §2.1| ) and explain how the network is formed ( §2.2| ). In § §2.3|2.4[ we collect necessary 
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information from the algebraic graph theory Q, which will be used for describing the role of the network 
topology in dynamical phenomena analyzed in this paper. In Section |3] we introduce two scenarios (A and 
B) leading to distinct firing patterns produced by the single cell models and by synchronized networks of 
these models. In the remainder of this section, we analyze the first of these scenarios illustrated in Fig. [T] 
(a,b). First, in LenrniajSTTJ we show that in the single cell model, bursting can be destroyed with small 
noise (Figure [T^). This counter-intuitive result relies on the presence of a slow timescale in the dynamics 
of bursting. Even very small noise can have significant qualitative effects on the dynamics of a slow-fast 
system, because during the (long) periods of slow evolution there is sufficient time for large deviations 
(which are extremely unlikely on time intervals of order 0(1)) to develop. The proof of the lemma uses 
large deviation estimates and is adapted from 1 17 |. Thus, given a deterministic model of a bursting cell, one 
can switch its dynamics into irregular spiking by adding noise. We then show that when many such cells 
are coupled together, the effects of noise weaken and the bursting of the underlying deterministic model 
reemerges (Figure [T])). 

The analysis of the coupled system proceeds in two steps. First, we use the center-manifold reduction 
lf37ll to approximate the coupled system near an excitable equilibrium by a simpler lower-dimensional system 
of ordinary differential equations. Second, thanks to the gradient structure of the reduced problem, we can 
accurately estimate expected time that a trajectory of the fast subsystem of the coupled system spends near 
the excitable equilibrium. We show that this time is much longer for the network model than for the single 
cell one. This analysis (based on large deviation estimates ifTSll ) yields one way of quantitative description 
of denoising. 

In Section |4j we present a complementary method based on a slow-fast decomposition. We show that 
when the coupling is strong, network dynamics near the excitable equilibrium splits into two modes: fast 
synchronization and ultra-slow noise-driven escape from the basin of attraction of the equilibrium along a 
low-dimensional synchronization subspace. The results of this section yield valuable insights into synchro- 
nization properties of the coupled system. In particular, we estimate the rate of convergence of trajectories 
to the synchronization subspace and the stability of the latter against random perturbations. The estimates 
show explicitly the contribution of the structural properties of the network to stability properties of the 
synchronization subspace. 

In Section [5j we take a look at denoising from a slightly different angle. Specifically, we study the 
linearization of the coupled system near an excitable equilibrium directly without invoking center manifold 
reduction. The analysis shows that when the dynamics of the individual cells lives in multidimensional phase 
space (unlike integrate-and-fire models or those of excitable cells near a saddle-node bifurcation), unless the 
coupling is full rank (see 114811 for the definition of full versus partial rank coupling) denoising should not 
be expected. These results show that the proximity to a saddle-node bifurcation in the fast subsystem of a 
square wave bursting neuron model, which makes the dynamics near the excitable equilibrium essentially 
one-dimensional, is critical for observing distinct dynamics generated by the single cell and network models. 
Therefore, the bifurcation structure of the bursting cell models like those used in (651 IIH and in this paper 
is important for the interplay between electrical coupling and noise in shaping the dynamics of the coupled 
network. Finally, the results of this study are summarized in Section [6j 
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Figure 2: The bifurcation diagram of the fast subsystem p.6|). Depending on the location of the null surface 



S, the full deterministic system (2.12)o and (2.13 ) can be in one of the three regimes: bursting (a), spiking 
(b), or excitable (quiescence) (c). The timeseries in (d) and (e) illustrate bursting and spiking respectively 



(see ^2.1 and Appendix B for details). 



2 The model 



In this section, we introduce a model that will be studied in the remainder of this paper. First, we formulate 
the assumptions on the single cell model and describe its main dynamical regimes. Next, we introduce the 
coupled system. At the end of this section, we review some notions and facts from the algebraic graph 
theory [2J, which will be useful for characterizing the contribution of the network topology to dynamics. 



2.1 The single cell model 

In numerical experiments throughout this paper, we are going to use a conductance based model of a pan- 
creatic cell due to Chay 0. The analysis in the following sections does not depend on any specific 
details of this model and more general assumptions will be formulated below. However, we believe that it is 
instructive to start from a concrete model to make the biological interpretation of the analysis that follows 
transparent. With this in mind, following [5J, we introduce the a system of differential equations modeling 
electrochemical dynamics in the /3— cells: 



CmV = -Iion{v,n,u) + aiwi, (2.1) 

riooiv) - n . 

n = — \-(J2W2, (2.2) 

t{v) 

u = e{Ica{v) -ku). (2.3) 



Here v, n, and u stand for the cell membrane potential, gating variable, and the concentration of calcium 
respectively. lioni^^ ^5 u), Ica{^), ^00 (^) and t{v) denote nonlinear functions, which are used for modeling 
ionic currents. Cm denotes membrane capacitance. The small parameter e > multiplying the right 
hand side of the third equation reflects the separation of the timescales of the calcium dynamics and the 
fast variables v and n. The right hand sides of the first two equations also contain independent copies of 
Gaussian white noise wi^2, which account for the deviations from the deterministic dynamics due to various 
fluctuations ^6l[7Tl. For further details of ( |2.1| )-( |2.2| ) including the values of parameters, we refer the reader 
to Appendix A and 0. 



To describe the structure of (2.1 )-(2.3 ), it is convenient to rewrite it in a more general form 



y = eg{x,y), 



(2.4) 
(2.5) 



where x = (xi,X2)^ := y := u, and H = diag(ai, (72). Let us first consider the deterministic 



model ( |2.4[ )o and ( |2.5[ ), where the zero subscript indicates that the stochastic perturbation is set to zero. 



0. The fast subsystem associated with ( 2.4 )o and ( 2.5 ) is obtained by setting e 
y as a parameter: 

X = f{x,y). 



in ( 2.5 ) and treating 
(2.6) 



Under the variation of y, the fast subsystem has the bifurcation structure as shown schematically in Fig. 
[2^. Specifically, 



(PO) There exists yhc ^ ^ such that for each y < y^c. Equation ( |2.6| ) has an exponentially stable limit 
cycle of period T{y): 

Ly = {x = 0(5, y): 0<s< T{y)}. (2.7) 

The family of the limit cycles, L = UyKy^^ forms a cylinder in R^, which terminates at a homo- 
clinic loop at y = yhc l|321 (Fig. [2k). 



(EQ) There is a branch of asymptotically stable equilibria of (2.6), E = Uyyy^n ^y^-^y ^ ^ ^ (?/)}' 



which terminates at a saddle-node bifurcation at = ysn < Vhc (Figure pa). 



(LS) For each y G M, the cj— limit set of almost all trajectories of (2.6) belongs to Ly [j Ey. 



Deterministic models of bursting are well understood (see, e.g.,|[2Tl|4Tl|5Tl|58]|). For small e > 0, (2.4)o 



and ( |2.5| ) features three main regimes: bursting, spiking, and quiescence (or excitable). In the former, the 
trajectory alternates between drifting along the cylinder L foliated by periodic orbits of the fast subsystem 
and the curve of equilibria (see Fig. [2^,d). Alternatively, (2.4)o and (2.5) may have a stable limit cycle 
near L (spiking) or a stable fixed point near E (excitable). The latter two regimes are illustrated in Fig. [2] 
(b,e) and (c) respectively. The analytical conditions for bursting, spiking, and excitable regimes are given in 
Appendix B. 
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2.2 The electrically coupled network 



Next, we consider a gap-junctionally coupled ensemble of n cells, whose dynamics is generated by ( |2.1| )- 
(2.3 ). In the coupled network. Cell z, z G [n] := {1, 2, . . . , n} receives current 

N 

I^^^9Y.a,,{v^^^-v^^), (2.8) 

from other cells in the network. Conductance aij > if Cell i and Cell j are connected and aij = 0, 
otherwise. Without loss of generality, we set an = 0, z G [n], and denote A = {cLij)- Nondimensional 
parameter ^ > is used to control the coupling strength. 



By including the coupling current p.8| ) into the models of individual cells ( |2.1| )-( |2.3| ), we obtain a dif- 
ferential equation model of the electrically coupled network 

Cmv^^ = -Iion(t'^'\n«,y«) +5 J^ai,-(w«) - + cjiw^'^^\ (2.9) 

y^"^ = e(/ca(i;»)-%»), z = l,2,...,n, (2.11) 
where VF^^^ = (it;^^'^), it;*^^'^))^ are independent copies of 2D Brownian motion. Using the notation, which 



we adopted for the single cell model in ( |2.4| ) and ( |2.5| ), we rewrite the coupled system in the following more 
general form: 

X = F{X,Y)-g{L®Ji)X + {In®Y.)W, (2.12) 
Y = G(X,y), (2.13) 



where 



and ® stands for the Kronecker product. 



2.3 The graph of the network 



Dynamics of the coupled system depends on the spectrum of matrix L appearing in the coupling operator 
in ( |2.12 ). The eigenvalues of L in turn depend on the structure of the graph of the network. Throughout 
this paper, we will repeatedly use the relation between the spectrum of L and the structural properties of 
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the network to study how the network topology affects its dynamics. To this end, we will need certain 
constructions and results from the algebraic graph theory, which we review below following 1(21. 



Let G = (V{G), E(G)) denote the graph of interactions between the cells in the network. Here, 
V{Q) = {^1, ^2, • • • , ^n} and E{Q) = {ei, 62, . . . , e^} denote the sets of vertices (cells) and edges (pairs 
of connected cells), respectively. Throughout this paper, we assume that G is an undirected connected graph. 

It is instructive to consider first the case when all nonzero conductances are equal to 1: 

1 Cell i and Cell 7 are connected, 

1 • (2.15) 
otherwise. 



As before, we set the diagonal elements of A to zero, an = 0. Matrix A = (a^j) in (2.15) is called the 
adjacency matrix of G and 

L = diag(deg(i;i), deg(i;2), . . . , deg{vn)) - A, (2.16) 

is called the Laplacian of G- By deg{vi) we denote the degree of Vi, i.e., the number of edges incident to Vi. 
Alternatively, the graph Laplacian can be defined by 

L = H^H, (2.17) 

where H G MJ^^^ is the coboundary matrix of ^ Q. The definition of H uses an orientation of the edges 
of G' For each edge ej = {vj^ , vj^) G V{G) x V{G), we specify the positive and negative ends; e.g., let vj^ 
be the negative end of ej^. Then the coboundary matrix of G is defined as follows (cf. Q) 

{1, Vj is a positive end of e^, 

-1, Vj is a negative end of ei, (2.18) 
0, otherwise. 

Definitions ( |2.16| ) and ( |2.17| ) are equivalent (cf. f2]). From either of these definitions, it is easy to see that 

Ai(L) = is an eigenvalue of L. If ^ is a connected graph then the zero eigenvalue is simple and all other 
eigenvalues are positive (cf. ll20ll ) 

= Ai(L) < A2(L) < . . . < Xn{L). (2.19) 

Remark 2.1. Following a common in the algebraic graph theory convention, we will refer to the eigenvalues 
of L as the eigenvalues of G- 

The second eigenvalue a(^) = A2(L) is called the algebraic connectivity of G, because it yields a 
lower bound for the edge and the vertex connectivity of G |20|. The algebraic connectivity is important 
for a variety of combinatorial, probabilistic, and dynamical aspects of the network analysis. In particular, 
it is used in the studies of the graph expansion ll33l , random walks |3 |, and synchronization of dynamical 
networks ||23l|47l[55l. Below, we show that a(^) determines the rate of convergence to synchrony in the 
coupled model ( [2l2| ) and ( |2l3 ). 

Another spectral function of G, which will be useful in the analysis of the coupled system is the total 
effective resistance of G (cf. ll40l[73]| ) 

n 

n{g)^nY,^]\L). (2.20) 

i=2 
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Figure 3: Symmetric (a) and random (b) degree— 4 graphs defined in Example 



2.3 



For electrical and graph-theoretic interpretations of TZ{Q) as well as for numerous applications, we refer the 
reader to |[3l[25]|. 

All definitions and constructions, which we reviewed above, naturally extend to cover nonhomogeneous 
networks with different conductances a^j. To this end, we define conductance matrix 

C = diag(ci,C2,...c^), (2.21) 

where q > is the conductance of edge Ci^i G [m]. The graph Laplacian of the weighted graph Q = 
{V,E,C) is defined by 

L{g) = H^CH. (2.22) 

The algebraic connectivity and the total effective resistance of Q are defined through the EVs of L{Q) as 
before. 



2.4 Examples of the network connectivity 

In this subsection, we present several examples of network connectivity, which will be useful for illuminating 
the role of the network topology in the dynamical phenomena analyzed in this paper. 

Example 2.2. One of the most common examples of local connectivity is a ID nearest neighbor coupling. 
It has been used in many studies of the coupled ensembles of /3— cells (see, e.g.,H65\ \64\l ). We use the ID 
nearest-neighbor coupling in most numerical experiments throughout this paper 

In this configuration, each cell in the interior of the array is coupled to two nearest neighbors. This 
leads to the following expression for the coupling current: 

lij) = - z;W) + g(v^^-^^ - v^^^), j = 2, 3, . . . , n - 1. 

The coupling currents for the cells on the boundary are given by 

/(I) ^^(^(2) _^(l)) =^(i;(^-l) -i;W). 
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The corresponding graph Laplacian is 
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(2.23) 
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The undirected graph corresponding to the ID nearest neighbor coupling scheme is called a path and is 
denoted by Pn 



Clearly, the total number of edges in the graph (connections in the network) is one of the important fac- 
tors shaping the dynamics of the coupled system. However, the way how these connections are distributed, 
i.e. the connectivity of the network, is also important as the analysis of the following example will show. 

Example 2.3. Consider two graphs on n nodes of degree d — A (i.e., each node in each of these graphs has 
precisely four edges incident to it). Such graphs are referred to as (n, d) —graphs. Thus, each of these graphs 
has 2n edges. To illuminate the role of connectivity, we assign to these graphs two different connectivity 
patterns as shown schematically in Fig. |j] The graph in Fig. ^ has symmetric connections. The edges of 
the graph in Fig. ^ were selected randomly. 

Specifically, the second graph was generated using so-called permutation model of a (n, d) random 
graph of even degree d = 2d In this model, one chooses at random d permutations 

7ri,7r2, . . . ,7rj 

in the symmetric group Sn- Then the edges between n vertices i^i, ^^2, • • • , ofQ are generated as follows 

E^{{vj,^^{vj)) : j G [n], i G [d]}. 



Spectral properties of the random graphs similar to the one constructed in the previous example have 
important implications for analyzing dynamical phenomena like synchronization and denoising in large 
networks. In the remainder of this subsection, we will discuss several facts about the spectra of random 
graphs that are particularly relevant to the analysis that follows. Our review is very brief and we refer an 
interested reader for more information and extensive bibliography to an excellent survey by Hoory et al |[33]| . 



First, we note that the random graph constructed in Example [23] is a (spectral) expander, which means 
that for some positive a 

)^2{Qn) > uniformly in n G N, (2.24) 



where n stands for the power of the set of vertices 1331 For those readers who have not seen this 
condition before, we note that this property fails to hold for any family of lattices (see Section 5 in BTll for 



a related discussion). In particular, for a path on n vertices (cf. Example 2.2), the second eigenvalue 



X2{Pn) = 4sin2 = 0(n-2) (2.25) 

tends to zero as n oc. Thus, the existence of a uniform bound for the second eigenvalue postulated by 



2.24 is a special, if not counter-intuitive, property. Nonetheless, it holds for random graphs. There are also 
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known explicit (nonrandom) algorithms generating expanders, including the celebrated Ramanujan graphs 

iain. 

In many applications, it is desirable to have a large bound on the expansion constant a. However, for 
(n, d) graphs, \2{Qn) can not exceed the Alon-Bopana bound g{d) = d — 2\/d — 1 031 . A remarkable 
property of the family of (n, d) random graphs is that for n » 1 they posses nearly optimal expansion 
constant with overwhelming probability. In particular, it is known that for any e > 

Prob \\2iQn) >d- 2Vd-l - e| = 1 - o^(l) Ve > 0, (2.26) 

where Qn stands for the family of (n, (i) random graphs of degree d>3 and n » 1 |[T9ll . 



3 Transitions to bursting in the coupled model 

Having reviewed bursting in deterministic systems and the definition of the coupled network, we now turn 
to the main theme of this work - the roles of noise and electrical coupling in shaping firing patterns of 
the coupled system ( |2.12| ) and ( |2.13| ). Fig. [T] shows transitions in the coupled system's dynamics observed 



over long intervals of time upon increasing the coupling strength. In the first case, which in the sequel we 
will refer to as Scenario A, irregular spiking patterns for weak coupling are transformed into fairly regular 
bursting patterns for sufficiently strong coupling (see Fig.[T^,b). In the second case. Scenario B, robust very 
irregular bursting becomes synchronous spiking when coupling strength is increased. 

We show that at the heart of both transitions lies denoising, the mechanism responsible for greatly 
diminishing the effects of noise on network dynamics. Below, we discuss both scenarios illustrating them 
with numerical results. For Scenario A, we also develop analytical estimates characterizing denoising. Our 
goal is to show how statistical features of the firing patterns depend on the coupling strength, excitability, 
and network topology. The analysis of Scenario B can be developed along the same lines, but it requires 
certain additional techniques for dealing with the analysis of trajectories near a limit cycle. These extensions 
will be considered in the future work. 



3.1 Scenario A: The single-cell model 



In this and in the following subsections, we discuss Scenario A in more detail. We start with the single-cell 



model (2.4) and d23 



Assume that the deterministic model (2.4)o and (2.5) is in the bursting regime, i.e., a typical trajectory 



alternates between drifting along a curve of stable equilibria of the fast subsystem, E, and a cylinder of 
limit cycles, L (see Fig. [2^). Let us focus now on a slow evolution alon g Since y = 0(e) (cf. (2.5)), 
on time intervals o(e~^) long, the dynamics of the slow-fast system (2.4)-(2.5) is approximated by the fast 
subsystem 

X = f{x,y) + S^, (3.1) 

where y G {ysn^Vhc) is fixed, x = (xi,X2)"^ G M^,S = diag(cri, cr2), cri,2 > 0, / : : 
smooth function, and w = {wi^W2) is a standard Brownian motion in M?. 



IS a 
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Figure 4: The phase portrait of the fast subsystem, (a) The basins of attractions of the stable fixed point 
Ey and limit cycle Ly are separated by the stable manifold of the saddle fixed point, (b) The trajectory of 
the randomly perturbed model \2.12) and ( |2.13| ) switches between the neighborhoods of the two attractors. 
Note that the transitions take place along the unstable manifold of the saddle. This shows that the center 
manifold reduction used in Sections |3] and |4] adequately describes the dynamics observed numerically. 



The frozen system ( |3.1| )o is bistable. It has two co-existing attractors: the stable equilibrium Ey and 
the limit cycle Ly (see Fig. [4] a). A trajectory of the randomly perturbed system (3.1 ) alternates between 
the basins of Ey and Ly, B{Ey) and B{Ly), separated by the stable manifold of the saddle point O (see 
Fig.|4|3). In fact, most of the time it spends in small neighborhoods of Ey and Ly. The plot in Fig.[4]3 shows 
that a typical trajectory leaves the neighborhood of Ey along the weak stable manifold of the sink, i.e., the 
dynamics near Ey is effectively one-dimensional. We will use this observation to simplify the analysis of 
the coupled system by reducing it (via the center manifold theorem 181) to a simpler system. 

The following exit problem is instrumental for understanding the effects of noise on the dynamics of the 



bistable system (3.1 ). For a trajectory of (3.1 ), which starts from a deterministic initial condition x(0) = 
xq G B{Ey), define the first exit time from B{Ey) by 



T{Ey,XQ) = inf{t > : x{t) ^ B{Ey)}. 



(3.2) 



Below we show that one can choose a(e) > such that with overwhelming probability for small e > 0, 

T{Ey,XQ) = o(e"^) 



and at the same time cr(e) as e 



can make the trajectory of the frozen system (3.1 ) leave B{Ey 
amount thus destroying bursting. 



(see Lemma 3.1). Therefore, with arbitrarily small noise one 
in time insufficient for y to change by 0(1) 



In Fig. [5^, we show a typical trajectory of ( |2.1[ )-( |2.3| ) superimposed on the bifurcation diagram of the 
fast subsystem. Note that the trajectory can not advance far enough along the curve of stable equilibria 
E, because it is forced to jump to the vicinity of L by noise. Thus, it lands on L near y^c every time and 
does not have enough room to generate many spikes before it reaches the right end of L and jumps back 
to E. This results in very irregular bursting patterns with very few spikes in one burst (see Fig. [5]3,c). In 
theory, for sufficiently small e > and small a > one can make the probability of clusters of 2 and 
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more spikes in the timeseries of ( |2.1[ )-( |23| ) arbitrarily low. This means, that one can transform bursting 
into irregular spiking by adding small noise. Showing this numerically, however, requires integrating stiff 
stochastic differential equations over very long time. Thus, we did not strive to achieve irregular spiking 
in our experiments resorting to irregular bursting patterns with very few spikes in Fig. |5]3, which already 
illustrate this effect. 

The noise-induced irregular firing pattern shown in Fig. ^ is characterized by the approximetelly ge- 
ometric distribution of spikes in one cluster (see Fig. [5]:). The geometric distribution has its origins in the 
exponential distribution of the first exit time r{Ey^ xq) II13I , which implies roughly that the distance from 
the landing point on L, yi, to the right hand end of L, y^c, Vi — Vhc is distributed approximately exponen- 
tially. The exponential distribution of yi — yhc translates into the geometric distribution of the number of 
spikes in one burst. Later we will see that this distribution is qualitatively different for the coupled system 
(see Fig. [5]0. The distinct probability distributions in Fig. [5]: and Fig. |5]P corresponding to different values 
of the coupling strength show that the transition in the network dynamics has taken place. 

In conclusion, we note that while the fact that we were able to destroy deterministic bursting with noise 
may not seem very surprising, the possibility of doing this with small noise is far from obvious. Note that 
the equilibria of the fast subsystem, Ey, for y near y^c are extremely stable. The slow-fast structure of the 
vector field is the key to this important effect. Extremely slow evolution along E gives the trajectory of the 
frozen system enough time to develop large deviations (necessary to leave B{Ey)), which are highly unlikely 
on time intervals 0(1) long. This is a general mechanism by which adding noise to slow-fast systems may 
create qualitatively new dynamical regimes ifTTll . 

Before we move on to discuss the coupled system, we prove the following lemma, which provides an 
estimate of the noise intensity sufficient for destroying bursting. 



Lemma 3.1. Let k > and S = diag (a, ka) be a matrix defining the noise intensities in {3.1). Then there 
exists Cl > such that for every C2 > Ci and 

a = (3.3) 



a trajectory of {3.1) with initial condition x(0) G B{Ey) leaves B{Ey) in time o(e ^) with probability 



converging to 1 as e ^ 



Proof. By rescaling X2 in ( |3.1| ), one can always achieve A: = 1, which we assume without loss of 
generality. 

Thanks to the large deviation theory, we have the following estimate for the first exit time T(£^y, xq) (cf. 
Theorem 4.4.2, HI): for any /i > 0, 

limP,, {exp{(K, - h)a-\e)} < T{E{y),xo) < exp{{Vy + h)a-\e)}} = 1, (3.4) 
where positive constant Vy is the minimum of the quasipotential associated with the randomly perturbed 



system (3.1 ). The definition and the properties of the quasipotential can be found in [ 18 |. 



Fix < a < 1 and take 

-,2_Vy + h 



c 



^ 1-a 
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Figure 5: Scenario A. The deterministic model ( |2.1| )-( |2.3| ) is tuned to the bursting regime (a,b). Weak noise 
can destroy bursting in the single cell model creating irregular spiking pattern (cf. Lemma [STT] ). The plots 
of trajectory of the coupled system superimposed on the bifurcation diagram (d) and the corresponding 
timeseries (e) show that bursting is recovered in the network thanks to denoising. The transition to bursting 
can be clearly seen from the normalized histograms of the number of spikes in one cluster or burst. The 
geometric distribution corresponding to (a,b) is transformed to the Gaussian distribution in (c,d). In these 
numerical experiments, we used system ( |2.1| )-( |2.3| ) with 50 oscillators coupled through the nearest neighbor 
coupling (see Example [Z2| ) with the coupling strength g = 5000 and other parameters specified in Appendix 
A. 



The combination of (3.4) and (|3.3|) implies 



T{E{y),xo) < e 



-Vy + h 



(3.5) 



with probability tending to 1 as e ^ 0, provided < h < Vy. 

□ 



3.2 Scenario A: The coupled system 

Next, we turn to the analysis of the coupled system. We want to understand how bursting is recovered for 
larger values of the coupling strength (see Fig. [5]d-f). Thus, we consider the coupled system ( |2.12| ) and 
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( |2.13| ). Below, we use two simplifying asuumptions, which let us avoid certain technical details, which are 



peripheral to the mechanism analyzed below. First, we focus on the fast subsystem ignoring 0(e) changes 
in the slow variables: 

X = Y) - g{L Ji)X + (/^ (3.6) 
where X = (x^^), x^^), . . . , x^""))"^ G M"^ x • • • x = R^^. Furthermore, it can be shown that with the 



diffusive coupling like in (3.6), y[s synchronize and remain close after some initial transients, provided the 
coupling is sufficiently strong (see ||48l|49|). Thus, we set the frozen slow variables for all subsystems to the 
same value y G (y^n, Vhc) so that 



For the coupled system (3.6), we consider the problem of exit from the basin of the rest state Ey = 
Ey X • • • X Ey e M?^. To this end, we define 

r(B(Ey),Xo) = inf{t > : X(t) ^ B{Ey)}. 



The key difference between the coupled system p.6| ) and the individual subsystems ( |3.1| ) is that the 



former is much less susceptible to the effects of noise. The robustness of the coupled system to noise 
manifests itself in the disparity of the exit times: 

T{B{Ey),Xo):^r{B{Ey),xo). 

More precisely, the asymptotic relation between the exit times corresponding to the single cell and coupled 
models is given in the following theorem. 



Theorem 3.2. Suppose that (3J_) is close to a nondegenerate saddle-node bifurcation^ Then for some 



(Jo > and go > 0, and for all < a < ao and g > go the following asymptotic relations hold: 

t(;B(£;^),xo) X exp and t(/3(E^), Xq) x exp j^^j , (3.7) 



where x denotes logarithmic asymptotics (cf jjjSA^ ), n is the number of the cells in the network, and C3 is 
a positive constant independent of a and n. 



Remark 3.3. Relations in ( |3.7| ) show that for the level of noise chosen in Scenario A, by taking sufficiently 
strong coupling with sufficiently many cells in the network, one can make t(;B(E^), xq) longer than the 
time necessary for y to reach the vicinity of ysn- This means that the level of noise, which prevents a single 
cell model from bursting, does not affect bursting in the coupled system (compare Fig. [5] a and d). 

Proof The proof of the theorem follows from the analysis of a closely related model in [|55il . For com- 
pleteness, we outline the main steps of the proof and refer the interested reader to |55j for further details. 
The proof consists of several steps, the main of which are the center manifold reduction of the single cell 



and coupled models near the saddle-node bifurcation (cf. (3.12)) and the variational interpretation of the 



reduced systems (cf. (3.15)) 



^The nondegeneracy conditions are specified in the proof of the theorem. 
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1. By our assumptions on the fast subsystem, the Jacobian i?/(0, 0) has a ID kernel. Denote 

e G kerL>/(0,0)/{0} and p G ker(L>/(0, 0))"^ such that p'^e = 1. (3.8) 
In addition to standard nondegeneracy and transversality conditions for a saddle-node bifurcation 

1 d'^ 

ai = ^g:^P^f(ue,0)\u=o y^O. (3.9) 
a2 = ^pV(0,m)I/.=o 7^0, (3.10) 

we assume that 

a^^p^Je^O, (3.11) 

where Ji is the matrix involved in the coupling operator (see p.6| )). Condition p.ll| ) guarantees that 
the projection of the coupling onto the center subspace is not trivial. Without loss of generality, we 
assume that nonzero coefficients ai^2,3 are positive. 

2. Under the conditions in 1., near the bifurcating equilibrium the coupled model can be reduced to an 
n— dimensional slow manifold. The reduced system (after appropriate rescaling of the dependent and 
independent variables and dropping higher order terms) has the following form: 

z = - In- jLz-i- aW, (3.12) 

where 7 is the rescaled coupling strength, ln = (l,l,...,l)^GM^, and W is the standard Brownian 
motion in R^Q 

Similarly, the single cell model is reduced to the following ID equation: 

C^C^-l + awt. (3.13) 



3. We recast ( |3.13| ) and ( |3.12| ) as randomly perturbed gradient systems. The former is rewritten as 

^ = -^'{^) + aw, $(C) = ^ + C-y. (3.14) 



For the latter, we use the structure of the coupling matrix ( |2.22| ), to reduce it to the following form 

• - _^ 

dz 



U^{z) + aW, U^{z) = ^{Lz, Lz) + ^(^^)' (^.15) 



1=1 



where L = \fCH and \/C = diag(y/ci, y/c2, . . . , ^/c^) is a square root of the nonnegative definite 
conductance matrix C (cf. |2.21| ). 

4. The large deviation estimates yield the logarithmic asymptotics for the exit times associated with the 
stable fixed point of ( [314] ) (cf. Theorems 2.1 & 3.1 in 



lim P {exp{(2$(l) - h)cj-^} < r{B{-l), Co) < exp{(2$(l) + h)a-^}} = 1, V/i > 0, 



(3.16) 



^ Throughout this paper, we use W to denote multidimensional Brownian motion in Euclidean spaces of different dimensions 
associated with the coupled systems. We reserve w to denote Brownian motions used to perturb single cell models. The dimension 
of the space for stochastic processes will be clear from the context and should not cause any confusion. 
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where $(1) = 4/3 is the value of the potential at the barrier at( = 1. Thus, 



t{B{-1),Co) X exp , 2 



(3.17) 



which shows the first relation in (3.7). 



5. The (deterministic) coupled system ( |3.12| )o has a stable fixed point at z = — In- The basin of attraction 
of this fixed point is bounded by 



dD^lj A, A = = (zi,Z2, ...,Zn)eR^: {3i e [n] : = 1) & (z^ < 1, G [n])} (3.18) 

i=l 

The estimate of the exit times from D follows from the analysis of the minima of the potential function 
Uj{') on dD. In ||55]| , we prove that for 7 > 2Ai(L^)~^, U^{z) achieves its minimal value on dD at 

Z In. 

4:71 

:= U^{ln) = y . (3.19) 

(Here, Ai(L^) denotes the smallest eigenvalue of matrix L^, obtained from L by deleting the first row 
and the first column.) Knowing the minimum of the potential function on the boundary of the basin of 
attraction of the stable fixed point, the large deviation estimates for the randomly perturbed gradient 
systems yield 



(3.20) 



□ 



3.3 Scenario B 

In this scenario, the deterministic single cell model is tuned to be in the spiking regime, i.e., it has a stable 
limit cycle Ly^ , yc < yhc in the vicinity of L (see Fig.l2b,e). When noise is added to the modeling equations 



(2.1 )-( |2.3[), it forces the trajectory to leave the basin of attraction of the limit cycle B{LyJ once in a while 
(see Fig. |6^), producing irregular bursting (Fig. |6]). This mechanism of irregular bursting was studied in 
detail in |31 1. In particular, it is shown in |31| that the number of spikes in one burst has asymptotically 
geometric distribution with parameter p ^ exp{C4 cr~^} for some positive constant C4. 

The normalized histogram in Fig. [6]: shows that with the level of noise chosen for this experiment, 
the system generates very long irregular bursts. In the second row in Fig [6jl-f, we show the results for 
two coupled cells and the coupling strength g = 50. Taking just two cells already changes the statistics 
of bursting significantly (compare Fig. |6]c and f). For three coupled cells, the bursts become even longer 
preserving the geometric distribution (see Fig. [6]i). For ten coupled cells and the coupling strength g = 500, 
we were unable to detect a single termination of spiking activity. While the coupled system is technically 
still in the bursting regime, the bursts become so long that for all practical purposes, one can speak about 
effective transition to spiking (see Fig.|6]g,h). 

The changes in statistical properties of the firing patterns generated by the coupled system for increasing 
values of the coupling strength are due to denoising. Here, electrical coupling diminishes the effects of noise 
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Figure 6: Scenario B. The underlying deterministic system is in the spiking regime (see Fig. |2}3). Adding 
small noise transforms spiking into irregular bursting shown in plots (a,b). Coupling two cells together al- 
ready increases the duration of bursts significantly (d,e). For ten coupled cells the duration of typical bursts 
is extremely long so that for all practical purposes (g,h), it can be considered as a transition to spiking. The 
normalized histograms for the number of spikes in one burst are plotted for a single cell, two, and three cou- 
pled cells in c, f, and i respectively. We were unable to numerically detect any bursts for the 10 coupled cells 
due to their extremely long duration. All histograms show approximately geometric distribution confirming 
the large deviations nature of the bursting patterns. The parameters of the geometric distribution (i.e., the 
expected value of the number spikes per burst) increases significantly as we go from one cell to three cell 
network (c,f,i). 
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on synchronous limit cycle oscillations and results in longer times that the trajectory of the coupled system 
spends in B(LyJ for larger values of g. Qualitatively, this is the situation that we have already analyzed 
for Scenario A. However, the study of the Scenario B involves certain additional technical details needed 
for extending the analysis to systems of coupled limit cycle oscillators (see [49,48]). We will consider the 
analysis of Scenario B in the future work. 



4 Synchronization 



In this section, we study synchronization in electrically coupled network ( |2.12| ) and ( |2.13| ). Our motivation 



is twofold. On the one hand, synchronization is an important aspect of the network dynamics, because 
for sufficiently strong coupling the network activity becomes synchronized. In particular, all numerical 
results shown in this paper for the coupled system feature synchronous activity. On the other hand, we show 
that the mechanism of synchronization is closely related to that of denoising. In particular, the estimates of 
stochastic stability of synchrony, which we derive below, reveal the contribution of the network topology and 
the strength of connections distribution to synchronization and denoising. We show that networks with larger 
algebraic connectivity and smaller total effective resistance are more effective for implementing denoising. 
To illustrate this point, we use random and symmetric degree-4 graphs defined in Example [23] Our analysis 
and numerics show that random graphs and, more generally, expanders have good synchronization and 
denoising properties. The analytical approach, which we develop in this section, complements that of the 
previous section. 



4.1 The new coordinates 

To study synchronization, we introduce a special system of coordinates, which to leading order decouples 
the two principal modes of the system's dynamics: fast synchronization and ultra slow excitation due to 
noise. Remarkably, the latter is captured by a scalar stochastic ordinary differential equation. The slow- 
fast decomposition is used to show that synchronization takes place for sufficiently strong coupling and to 
quantify various aspects of the synchronized dynamics. 

The new coordinate system takes into account the structure of the network. To this end, we will need 
a spanning tree of ^, ^ = (^i^!), i.e., a subgraph of G with n = |^(^)| vertices, n — 1 edges, and 
containing no loops [2J. Having chosen the spanning tree let H E M(^~^)><^ denote the coboundary 
matrix corresponding to G- Matrix which we introduce in the next lenmia will be useful for constructing 
the new coordinates. 

Lemma 4.1. Let 

S = {HH^)-^/^H. (4.1) 

Then 

SS^ = In-i and 51n = 0. (4.2) 



Proof. Properties ( |4.2| ) follow from the definition ( |4.1[ ). 

□ 



19 



The synchronization subspace spanned by In coincides with the kernel of S, while the columns of 
form an orthonormal basis of In^. These two subspaces are important for studying synchronization, which 
motivates the following coordinate transformation 

W3 T]) G R""-^ X R, (4.3) 

where 

^ = Sz and r] = n-^ln^z. (4.4) 



Note that |^| = \Pi^±z\ measures the distance of the solution of the coupled system (3.12) to the synchro 



nization subspace corresponding to ^ = 0. Here, Pi^± stands for the orthogonal projector onto In • 



4.2 The slow-fast system 



Throughout this section, we work with the reduced (rescaled) system p.l2[ ), which we derived using the 
center manifold approximation of p.6| ) near the excitable equilibrium. In Section [3j we analyzed ( |3.12| ) 
using its gradient structure and the large deviations estimates. This time, we use a complementary approach 



by first identifying and then exploiting the slow-fast structure of ( |3.12| ). 
By projecting (3.12) onto In^ and In, we obtain 



v 



(-7L + 2r?I„_i)e + aSW + 0{\^\' 



n 



(4.5) 
(4.6) 



where matrix L = SLS^ is the unique solution of the matrix equation 

SL = LS. 

He re, W and w denote the Gaussian white noise processes in and respectively. In the derivation 
of (4.6k we use t he fact that In^Vi^ and n~^/'^w are identically distributed. For details of the derivation 
of ( 451 ) and (4^), we refer the interested reader to Lemma 5.3 of ll55l . About the reduced matrix L, the 
following is known: L is a positive definite matrix, whose spectrum consists of the nonzero eigenvalues of 
L (cf. Lemma 2.6, L42J)- Moreover, S maps the generalized eigenspaces of the nonzero eigenvalues of L 
bijectively onto those of L BTll . 



Equation ( |4.5| ) captures the dynamics along the orthogonal complement of the synchronization subspace 
In. Thus, it describes the process of synchronization. On the other hand. Equation ( |4.6| ) tracks the motion 
along the synchronization subspace. Since L is positive definite, in the strong coupling regime (7 » 1), it 



follows from (4.5) that the trajectory of the full system (4.5) and (4.6) relaxes to an 0{a) neighborhood of 
the synchronization subspace. In, at 0(7) rate. Equation (4.6) to leading order (with ^ p:^ 0) is a standard 
model of a particle in a potential well forced by noise. On the time intervals not exceeding the Kramer's 
time 0(exp{C5cr~^}) for some C5 > 0, 77 remains close to —1 for the most of the time. This sums up the 



qualitative dynamics of (|4.5|) and (4.6). In the remainder of this section, we study it in more detail. 
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4.3 The fast subsystem: synchronization 



The analysis of the fast subsystem elucidates several important aspects of synchronization in electrically 
coupled networks. In particular, we estimate the rate of synchronization in terms of the network connectivity. 
We then study robustness of synchrony to noise. 



The stability of the synchronization subspace is determined by the linear part of ( |4.5| ): 

^ = + aSW. (4.7) 

Assuming for simplicity the deterministic initial condition ^(0) = Co ^ R^~^, we compute the mean vector 
and the covariance matrix of ^ (t): 

EC(t) = exp{-7tL}Co, (4.8) 

cov i{t) = crM exp{-27L(t - u)}du = ^L~^ [in-i - exp{-2t7L} ) . (4.9) 
Jo 27 V / 



It follows from ( |4.8| ) and the geometric interpretation of C, that the rate of convergence to the synchro- 
nization subspace is set by the product of the strength of coupling 7 and the algebraic connectivity a(^): 



lEPi ^z{t)\ < C6exp{-7a(e)t}, 



(4.10) 



where P^^^ ' stands for the orthogonal projection onto In^ and Cq is a positive constant independent from 
7 and a. Therefore, networks with larger algebraic connectivity synchronize faster. There are many fine 
results in the spectral graph theory relating algebraic connectivity to various structural properties of the 
network (see [33 ] and references therein). These results can be used via ( |4.10| ) to elucidate the contribution 
of the network topology to synchronization properties of the coupled system. In particular, ( |4.10| ) shows 
that networks on random graphs and on expanders in general (cf. Example [23] ) admit a lower bound on the 
synchronization rate, which is uniform in the size of the network n. 

In the presence of noise, it is important to know how well the synchrony can withstand stochastic per- 
turbations. This leads to the question of stochastic stability. There are many ways for measuring stability 
of synchrony to random perturbations 130]. In this paper, we use the mean square stability, which pro- 
vides a natural metric for the problem at hand. Specifically, we are interested in transverse stability of the 



synchronization subspace. From (4.10), we know that 

^P^^±z{t) ^0, as t ^ oc. 
Thus, we next look at the second moments, i.e., we estimate the dispersion of the trajectories around In: 



n-l 



var P^^±z{t) = ^ var ^i{t) =: var ^{t). 



(4.11) 



From (4.9), we find that 



var ^{t) = Tr cov ^{t) = — Tr 

27 



(/n-i - exp{-2t7L} 
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Since L is positive definite, var ^ [t] has a finite asymptotic value 

_2 



var ^ := lim var^(t) 



27- 



-TrL" 



2n7 



(4.12) 



where TZ{Q) is the total effective resistance of the weighted graph of the network Q (cf. 2.20). Thus, vaT^ 
provides a convenient measure of stochastic stability of the synchronization subspace. The smaller the value 
of vaT^, the more stable the synchrony is. In the next section, we will use the asymptotic value of the 
variance to estimate the effectiveness of denoising. 



Estimate \A.\2) shows that the stochastic stability of the synchronization subspace is fully determined 
by the strength of coupling and the total effective resistance of the network. Similarly to the algebraic con- 
nectivity of the graph, the value of the total effective resistance TZ{Q) can be related to the structure of the 
graph and the weight distribution |25|. However, while the rate of convergence to the synchronization sub- 
space depends only on the leading nonzero eigenvalue of the graph Laplacian, the description of stochastic 
stability requires the entire spectrum of Q. The information about higher eigenvalues of Q in general is hard 
to obtain. However, as one can see from the following crude estimate of TZ{Q) 

.2 



n 



(4.13) 



graphs with larg er alg ebraic connectivity enjoy better bounds on the total effective resistance. In fact, for 
expanders, from (4.13 ) one gets a bound on TZ{Q) = 0{n^), which can not be improved (up to a multiplica- 
tive constant). Therefore, we expect that graphs with good bounds on algebraic connectivity, like a random 



graph in Example |2.3[ are robust against random perturbations. Below, we will illustrate this point with 
numerical results. 

Above we have used the analysis of the fast subsystem to gain useful information about synchronization 
properties of the coupled network. Further, with these results at hand one can easily understand the ID slow 



equation (4.6) and get a complete description of the slow-fast system (|4.5|) and (|4.6|). This is our next step. 



Since on time intervals not exceeding the Kramer's time, with high probability ^(t) = 0{a), we approx- 
imate (|4.6|) by 



V = fiv) + 



W. 



n 



(4.14) 



From ( |4.14| ), we estimate the time that a typical trajectory spends in the basin of the rest state 77 = —1 

8n 1 



t(-1,77o) - exp 



3a2J 



(4.15) 



Comparison of p.20| ) and ( |4.15| ) shows that the results of the fast-slow analysis are consistent with those of 
the previous section. Furthermore, through the estimates of stability of the synchronization subspace ( |4.10[ ) 
and ( |4.12| ) the analysis of this section elucidates the contribution of the structural properties of the network 
to its synchronization properties. 



5 The mechanism of denoising 



In this section, we continue to study denoising, the mechanism responsible for variability of activity patterns 
in electrically coupled network ( |2.9| )-( |2.10| ). The analysis of synchronization and denoising in the previous 
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sections relied on the proximity of the models of individual cells to the saddle-node bifurcation. To clarify 
to what extent the mechanism of denoising depends on this assumption and to elucidate the scope of its 
applicability, in this section, we analyze the coupled system around the stable fixed point without assuming 
anything about the bifurcation structure of the problem. 

Because the dependence of the single cell models on the slow variable is not essential for the mechanism 
of denoising per se, we omit to indicate it explicitly to simplify notation. Thus, in this section, we consider 

X = f{x) + ^w, xgM"^, (5.1) 

where / : ^ is a smooth function, it; is a standard Brownian motion in R^, and H G R^^^ is a 
nonsingular matrix. We assume that (|5.1|)o has a stable equilibrium at the origin, i.e.. 



/(0) = 0, -A^DfiO), (5.2) 
where := 0.5(A + A'^) > is a positive definite matrix. The coupled system has the following form 

X = F{X) - g{L J)X + {In T.)W, (5.3) 

where as before X = (xW, x^^), . . . , x(^))^ G R^ x R^ ^ R^^, F(X) = (/(xW), /(x^^)), . . . , /(xH))^, 
L — H^CH, and J is d x d matrix, such that is positive definite. Finally, W stands for the Brownian 
motion in R^^ The coupled system has an equilibrium at the origin in R^^. The linearization of (53 ) about 
this equilibrium is given by 

X = -gBX + {In ^)W, B^L®J + 5In®A,5^g-^. (5.4) 



In the previous sections, to quantify the effect of denoising we used the times of the first exit from the 



basins of stable equilibria of the single cell and coupled models (cf. ( |3.20| ) and ( |4.15| )). Since in this section 
we do not assume that the local systems are located near a saddle-node bifurcation, we no longer can rely 
on explicit estimation of the first exit times. Thus, we seek other means for measuring the effectiveness 
of denoising. To this end, we recall that in the previous section we found that the asymptotic value of the 
variance of the trajectories near the synchronization subspace to provide a convenient measure for estimating 
stochastic stability. We thus adapt it to our current purpose. Specifically, we use 

d 

maxvaFx^^^ = max lim var x^-\t) 

to measure the variability of the trajectories of the coupled system X{t) = {x^^\t) , x^'^\t) , . . . , x^'^\t)). 
We estimate the effect of denoising by comparing max^^[^] vaTx^^^ to vaTx, where x{t) is the solution of 



the local subsystem (5.1 ). We identify conditions which guarantee that 



maxvar x^^^ <C var x, 



and show how the former quantity depends on the coupling strength, network size and topology, as well as 
on the intrinsic properties of the local system ( |5.1| ). 

The following lemma is the key for understanding the mechanism of denoising. 
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Lemma 5.1. 

(i) ^ cr'^ f ^ n-1 \ 

maxvar x^'^ < — , , , , H , , , (5.5) 



where a = stands for the Frobenius norm ofT,, B is defined in \5.4\ , and Ai(-) denotes the smallest 
eigenvalue of a symmetric matrix. 

In the process of proving Lemma |5.1[ we derive the following estimates, which are of independent 
interest. 



Corollary 5.2. 



< —( + ? , \ ^ , (5.6) 
\EP^^^X{t)\ < Cexp{-^Ai(S^)t} and Y^P^^^X{t) ^ ^2g\^{B^y ^^'^^ 



Proof. The proof of Lemma 5.1 consists of the following steps. 



1 . First, we separate the dynamics along the synchronization subspace from that along its orthogonal 
complement. To this end, we switch to new coordinates 



X^{Y,Z),Y ^{S® Id)X and Z = {n-^lrJ O Id)X. (5.8) 



By multiplying both sides of (5.4) hy S ® Id and (n ^In® IdY ^ we obtain the equations for Y and 
Z 

Y = -gBY + {S®^)W, B^L(®J + 6{In-i(®A), (5.9) 
Z = -AZ + n-'^/'^T,w. (5.10) 

2. Recall that A'^ is a positive definite matrix and the smallest EV of A'^ denoted by Ai(^*). For the 
mean vector of the solution ( |5.10| ) subject to the deterministic initial condition Zq 



|E Z| = I exp{-tA}Zo\ < Cexp{-Ai(^^)t}|Zo|, (5.11) 
for some positive constant C. Next, we estimate 

var Z = n"^Tr (^J^ exp{(t - s)A}ET,'^ exp{(f - s)A'^}ds^ 

< n-^Tr(EE'^) / exp{2Ai(^")s}ds ^ (2n|Ai(^")|)-V2 as t ^ oo. (5.12) 
Jo 



3. Similarly, 



|IEr| = \exp{-tgB}Yo\<Cexp{-gXi{B'')t}\Yo\md (5.13) 
vary = \S^nl ^ {n-l)a^ 
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4. Using (5.8 ), we express X in terms of Y and Z 

X = (5^ Id)Y + (In Id)Z =: MiY + 
The definition of S implies 

mJMi = and M2^M2 = n/^. 

By ( |5.15| ) and ( |5.16| ), we have 

TrEXX^ = Tr{MiE(yy^)M7 + M2E(yy^)M2^} 
= Tt{mJMiE (YY^) + (YY^)} 
= Tr{E(yy^) + nE(ZZ^)}. 

Using the elementary properties of cov , ( |5.11| ), ( |5.13| ), and ( |5.17| ), we have 

vaTX — vary + nvaTZ. 

The combination of ^J^, (jSTT?]), and ([STTS]), shows ( [5^ . 

5. Since X = (x^^), x^^), . . . , x^), from (|5.15|) we have 



By noting 

from ( |5.19| ) we have 
Thus, 



x» = z + (RoWi(5^) /rf)y =: z + TV^y. 
A^TtV- = and = 0, 
Tr E x^xW^ = Tr E {zz'^ + yy^}. 
vaTx^^^ — vaTz + vary 



(5.15) 



(5.16) 



(5.17) 



(5.18) 



(5.19) 



(5.20) 



Using (5.20), (5.12), and (5.14), we obtain (5.5) 



□ 

We want to understand how the variability of the coupled system can be smaller than that of a local 
subsystem. Estimate ( |5.5| ) suggests a possible scenario. Note that the first term on the right hand side of 
( |5.5| ) can be made arbitrarily small by increasing n. The second term can be controlled by g provided that 
Ai(S^) remains 0(1) as ^ ^ 00. Thus, we need to understand the behavior of Xi{B^) for increasing values 
of ^. In the following lemma, we show that the latter depends on the coupling architecture and identify two 
cases of full and partial coupling, which are important in the context of denoising. 

Lemma 5.3. Let 

B ^L®J^5{In-i®A), 
where > 0, > 0, and {) < 5 <^ 1 (cf. {5.4\). If is a full rank matrix then 



\i{B') = \i{L)\r{r) + 0{5). 



(5.21) 
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Otherwise, let k — dimker > 0, choose {gi, • • • , ^fc} orthonormal basis for ker J^, and define an 
k X k matrix 



Then 



Xi{B')^5Xi{G) + 0{5^). 



(5.22) 
(5.23) 



Proofi For small 5 > 0, the EVs of B pe rturb smoothly from the EVs of L ® (cf. El). If > 
then Ai(L ® J^) = Ai(L)Ai(J^) and ( 5.21 ) follows. Otherwise, is an EV of L of multiplicity 
k{n — 1). We construct an orthonormal basis for the 0— eigenspace of L® 

(z2 - + n, n G [n - 1],Z2 ^ [fc], 



(5.24) 



where = (5|, ^l, . . . , and 5^ stands for the Kronecker delta. By the perturbation results for the 

multiple eigenvalues of symmetric matrices (cf. Appendix, O), we have 

Ai(S^) = 5Ai(G) + 0(52), 

where G is an (n — 1)0? x (n — l)d matrix 

(G)^^- = h]{In-i®A)hi. 

It is easy to see that G = In-i ® G with G defined in ( |5.22[ ). Therefore, the EVs of G and G coincide. This 
shows ( |5^ . 

□ 

In conclusion, we discuss the implications of the analysis of this section. There are two effects that 
contribute to denoising: averaging and synchronization. The former can be interpreted as a manifestation of 
the law of large numbers: the combined effect of independent stochastic processes acting on individual cells 
vanishes as the size of the network goes to infinity. It is captured by the first term on the right hand side of 
\5.5) . From \5.5) it follows that averaging depends on the network size, n, and the dissipativity of the local 
dynamics through \i{A^) but is independent from the properties of the coupling operator. The second term 
on the right hand side of ( |5.5| ), which captures synchronization, depends on the dissipativity of the coupling 
operator through Ai(S^), the coupling strength g, and the size of the network. Denoising takes place when 
the variance of the trajectories of the coupled system can be controlled by n and g. In fact, we can always 
make the first term in \5.5) arbitrarily small by taking n large enough. To control the second term on the 
right hand side of (5.5), we can use g provided lim^^oo > 0. Thus, one can make the right hand 

side of (5.5 ) arbitrarily small by increasing n and g provided Ai {B^) is bounded away from as ^ ^ 00. 



Lemma [53] identifies two cases important in the context of denoising and synchronization. Depending 
on the rank of J^, the leading eigenvalue of B^, Ai (B^), is either 0(1) or 0(5) for ^ 1. If the coupling is 
full rank then Ai(S^) = 0(1) and (5.5) shows that the variability of large coupled systems can be effectively 
controlled by the coupling strength and the network size. However, if the coupling is partial (rank < d), 
as in the model of gap-junctionally coupled network ( |2.9| )-( |2.lT] ), ( |5.5| ) becomes 

^2 



var x^^^ < 



max var x 

ie[n 



a 
~2 



1 



+ 



n 



1 



n\i{A') Ai(G^) + 0(5)y - 2Ai(A 



< 



+ 0{G'5n) 



(5.25) 



The second inequality in ( 5.25 1 follows from \522\ and the variational properties of the EVs of symmetric matrices 



Ai(G') — min 

a^Gker , \x\- 



x^A^x > min x^ A^x — \i{A^ 
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This shows that in the partial coupling case the mechanism of denoising may fail. Our numerical experiments 
and the analysis of ( |2.9| )-( |2.lT] ) show that denoising is at work despite the fact that the coupling is partial 
in this model. We are able to observe denosing in ( |2.9| )-( |2.lT] ) because of the additional structure of this 
problem - the proximity of the models of individual cells to the saddle-node bifurcation. This additional 
feature of the problem affords center manifold reduction, and the coupling of the reduced system is already 
full rank. This reconciles the results of this section with the previous analysis and our numerical results. 
This also highlights the importance of the bifurcation structure of the fast subsystem of the bursting cell 
model (2.1)-(2.3) for the mechanism of denoising. More generally, this discussion implies that while the 
mechanism of denoising can be already studied at the level of the linearized system, the structure of the 
nonlinear system is nonetheless important for the realization of denosing in concrete (biophysical) models. 



6 Discussion 



In this work, we have analyzed synchronous regimes in electrically coupled networks of bursting capable 
excitable cells. The individual cells in the network can be tuned to one of the three main activity states: pe- 
riodic spiking, bursting, or quiescence. Strong electrical coupling synchronizes activity across the network. 
However, in the presence of noise, the synchronous patterns generated by the network can differ qualita- 
tively from the patterns that the cells comprising the network exhibit in isolation. We have identified two 
scenarios of such behavior: Scenario A and Scenario B. In the former, the network of irregularly spiking 
cells generates very regular bursting provided the coupling is sufficiently strong. In the latter, the network 
formed from irregularly bursting cells is effectively switched to spiking once the coupling strength and the 
size of the network become sufficiently large. 

In constructing these scenarios featuring the disparity of the firing patterns generated by the single cell 
and coupled models, we used the large deviations type mechanisms to generate irregular patterns in the sin- 
gle cell models and denoising for shaping the patterns of collective behavior in the network. In particular, for 
irregular spiking in Scenario A, we used noise to perturb the system from slowly evolving stable equilibrium. 
A closely related mechanism was studied in the context of emerging regular dynamics in randomly perturbed 
slow-fast systems [ fTTl and self-induced stochastic resonance ||57l[T4l. The irregular bursting in Scenario 
B was organized by perturbing a system with a stable limit cycle via the mechanism proposed in |[3T1l . In 
both cases, the slow-fast structure of the single cell models was essential. The emerging firing patterns in 
both scenarios are very sensitive to the variations in the intensity of noise. Utilizing the ability of electrical 
coupling to synchronize the activity and to reduce the effects of noise on network dynamics EOlllllllll, for 
the coupled system we achieved effective control of the firing patterns by varying the strength of coupling 
and the size of the network. 

In this paper, we analyzed in detail the mechanism of denoising for coupled systems with (slowly evolv- 
ing) stable equilibria (Scenario A). We expect that the analysis of Scenario B can be done along the same 
lines by combining the techniques of Section [5] and local analysis near the synchronous limit cycle of the 
coupled system as in EH 133. We will address this problem in the future work. For Scenario A, we con- 
sidered two cases: one - when the individual subsystems are close the saddle-node bifurcation as in the net- 
work of model bursting neurons ( |2.12| ) and ( |2.13| ), our motivating example, and another without assuming 



the proximity to the saddle-node bifurcation. For systems near the bifurcation, we adapted two comple- 
mentary approaches from L55J . The former of which utilizes the gradient structure of the reduced system 
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Figure 7: Network connectivity and denosing. The trajectories of the coupled systems of 200 bursting cell 



models on a synmietric (a) and random (b) degree— 4 graphs respectively (see Example [23] ). The trajectories 
in (a) fill a larger area as in (b), which suggests that denoising is more effective in a network on a random 
graph. The corresponding timeseries are shown in (c) and (d); and the histograms for the number of spikes 
in one burst - in (e) and (f). The latter shows that the number of spikes in one burst generated by the random 
network has approximately Gaussian distribution tightly localized around 15. This is in a stark contrast with 
a much broader distribution corresponding to the symmetric network in (e). 
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to yield accurate estimates of the first exit times from the basins of the stable equilibria of the single cell 
and coupled models; while the latter relies on the algebraic graph theory techniques [47], which elucidate 
the contribution of the network topology to synchronization and denoising. Furthermore, we extended the 
analysis of denoising in systems of coupled integrate and fire models in |[50l to conductance based models 
with multidimensional phase spaces. We have shown that in this setting the coupling architecture is very 
important for implementing denoising. In particular, for the common in applications partial coupling case, 
we have shown that denosing generically does not take place, in contrast to what one might expect from the 
analysis of the coupled one-dimensional systems |50|. This result highlights the significance of the bifur- 
cation structure of the biophysical models of square wave bursting neurons for the realization of denoising. 
The proximity to the saddle node bifurcation makes the dynamics of the individual subsystems effectively 
one-dimensional and the coupling effectively - full rank, thus circumventing the problems for implementing 
denosing in partially coupled systems identified in Section [5] Therefore, the results of this study extend 
and complement the existing results characterizing denoising in electrically coupled networks |[50l [68l and 
highlight the importance of the bifurcation structure of the local dynamical systems comprising the network 
for implementing denoising. 

In analyzing the coupled system, we paid special attention to the role of the network topology in shaping 
the network dynamics. We have found that two spectral functions: the algebraic connectivity and the total 
effective resistance, feature prominently in the quantitative descriptions of synchronization and denoising in 
electrically coupled networks. The algebraic connectivity of the weighted graph of the network sets the rate 
of convergence to the synchronization subspace (cf. ( |4.10| )), while the total effective resistance is involved in 



the estimates of stochastic stability of the synchronous regime (cf. ( |4.12| )). These analytical estimates allow 
one to use many known results relating the spectra of the graphs and their structural properties (see ||9l |33l 
and references therein) to elucidate the contribution of the network structure to its dynamics. As an example 
of such application, we used the spectral properties of random graphs |[33l to show that networks on random 
graphs feature fast synchronization and robustness of synchrony to noise. They are also more effective 
for implementing denoising than their symmetric counterparts (see Example |2.3| ), as shown in numerical 
experiments in Fig. |7j 

Understanding dynamical mechanisms for different patterns of electrical activity in excitable cells and 
transitions between them has been long recognized as a fundamental problem in mathematical neuroscience. 
The results of this paper describe the transition from irregular spiking to nearly periodic bursting in networks 
of electrically coupled cells in the presence of noise. This transition was used in [|65l |64l to explain why 
pancreatic /3— cells burst in electrically coupled islets of Langerhans but not in isolation. Our results support 
and extend the previous analysis in |[65l|64l- We believe that the results and techniques of this work will be 
useful for understanding dynamics in many other biophysical models of gap-junctionally coupled networks. 

Acknowledgments. This work was partly supported by the NSF Award DMS 1 109367 (to GM). 



Appendix A. The parameter values used in the conductance-based model of 
the beta cell 

In the numerical experiments for this paper, we use a conductance based model of a pancreatic /3— cell due 
to Chay |5 |. Below, we collect the expressions of the nonlinear functions and the parameter values used in 
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The first term on the right hand side of ( |2.1| ) models the combined effect of sodium and calcium currents, 
iNa+Ca, the calcium-dependent potassium current, IxCa, delayed rectifier Ik and a small leak current, // 



lion = iNa+Ca + IxCa + Ir + h-. 



(A.l) 



where 



iNa+Ca 

Irc 

Ik 
k 



9Na+Camoo{vfhoo{v){Ei - v) , 

gi{Ei - v). 



The steady state functions used to model the ionic currents above are given by 



foc{v) 



af{v) + ^f{v) ' 



/ G {m,/i,n}, 



where 



= {).l{v + 25)(1 - exp{-0.1(i; + 25)})"\ = 4exp{-(i; + 50)/18}, an = 0.07exp{-0.05(i; + 50)}, 
= (1 + exp{-0.1(i; + 20)})"\ an = 0.01(i; + 20)(1 - exp{-0.1(i; + 20)})"\ = 0.125 exp{-(i; + 30)/80}. 

The time constant of the delayed rectifier is given by 

= (230K + /3,))-\ 

The values of the remaining parameters are summarized in the following table. 

Table A.l 



QNa+Ca 


ISOOs-i 


ENa+Ca 


lOOmy 


9k 


ITOOs-i 


Ek 


-15mV 


kc 




m 




El 


-40mV 


9KC 


12s-i 


Eca 


lOOmV 


e 


o.osmy-^s-i 


Cm 


IfiF/cm? 


a 





Appendix B. Dynamical regimes of the conductance-based model 



The geometric theory for singularly perturbed systems implies the existence of the exponentially stable 
locally invariant manifolds and L^, which are 0(e) close to E ^{{x^ y) : y > ysn + ^} and L ^{{x^ y) : 
y < Vhc ~ respectively, for arbitrary fixed S > and sufficiently small e > ||Tl|22l. Manifolds E^ 



and Le are called slow manifolds. For small e > 0, the dynamics of ( |2.4| ) and ( |2.5| ) on the slow manifolds is 
approximated by 



E, : 



y = eG{y), y<yhc-S, 



(B.l) 
(B.2) 
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where 



1 fT{y) 



{(p{s,y),y) ds 



(B.3) 



and T{y) stands for the period of the Hmit cycle Ly. De pendi ng o n the location of the null surface S = 
{(x, ?/) G : g{x^ y) = 0}, the slow-fast system (2.4) and (2.5) can be in one of the following states: 
bursting, spiking, and quiescent (see Fig.|2]). 

The following conditions on the slow subsystem yield bursting. For some c > independent of e, 



(SE) 
(SB) 



9{^{y)^y) < -c for y > ys 



G{y) > c for y < yhc- 



(B.4) 
(B.5) 



Under these assumptions, for sufficiently small e > a typical trajectory of ( |2.4| ) and ( |2.5| ) consists of the 
alternating segments closely following and and fast transitions between them (see Fig. [2^). The 
corresponding timeseries is shown in Fig.[2]l. 



Substituting (SB) with the condition that follows will switch (2.4) and (2.5 )o to a spiking regime. 



(SS) G(y) has a unique simple zero at y = ^ (Vsn, yhc)'- 

G{yc) = and G\yc)<0. 



(B.6) 



In this case, the Pontryagin-Rodygin theorem |56| yields the existence of an exponentially stable limit cycle 
Le(^c) of period T{yc) + 0(e) lying in an (9(e) neighborhood of Ly^, provided (SS) holds and e > is 
sufficiently small (see Fig.l2b,e). 



Finally, the slow-fast system (2.4)o and (2.5 ) is said to be in the excitable regime if it has a stable fixed 
point lying on (see Fig.l2b): 



(Q) 9{y) 9{^{y)^y) has a unique simple zero ^ty^y^e {ysn, yhc)- 

g{yq) = and g\yq) < 0. 



(B.7) 
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